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Abstract 

A macroscopic model to describe the dynamics of ion transport in ion channels is the 
Poisson-Nernst-Planck(PNP) equations. In this paper, we develop a finite-difference 
method for solving PNP equations, which is second-order accurate in both space and 
time. We use the physical parameters specifically suited toward the modelling of ion 
channels. We present a simple iterative scheme to solve the system of nonlinear equa- 
tions resulting from discretizing the equations implicitly in time, which is demonstrated 
to converge in a few iterations. We place emphasis on ensuring numerical methods to 
have the same physical properties that the PNP equations themselves also possess, 
namely conservation of total ions and correct rates of energy dissipation. We describe 
in detail an approach to derive a finite-difference method that preserves the total con- 
centration of ions exactly in time. Further, we illustrate that, using realistic values 
of the physical parameters, the conservation property is critical in obtaining correct 
numerical solutions over long time scales. 
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1 Introduction 



The Poisson-Nernst-Planck(PNP) equations describe the diffusion of ions under the effect of 
an electric field that is itself caused by those same ions. The system is created by coupling the 
Nernst-Planck equation (which describes the diffusion of ions under the effect of an electric 
potential) with the Poisson equation (which relates charge density with electric potential). 
This system of equations has found much use in the modelling of semiconductors. [H] Al- 
though the Poisson-Nernst-Planck equations were applied to model membrane transport for 
longer than they have been employed to model semi conductors [17], the use of the system 
to model the behavior of the internal mechanics of these transport processes is much more 
recent. [5] 

The system of PNP equations and its related models have been the subject of much 
study and numerical simulation. A recent advancement in this field was the application 
of energy variational analysis and density functional theory to modify the PNP system to 
accommodate various phenomena exhibited by biological ion channels. See p2] and the 
references therein. 

The computer simulations of the Poisson-Nernst-Planck models are able to capture the 
transient, dynamical behavior of the system, and the numerical schemes employed are quite 
varied. Cagni et al. (2007) [2] discretized the PNP in two dimensions using a second-order ac- 
curate finite difference method with central differencing in space and Crank-Nicolson scheme 
in time, and simulated an ion channel subjected to time-dependent perturbations. Nan- 
ninga (2008) [15] studied a nerve impulse using a similar finite difference scheme as in [2] 
but in three dimensions, notable in that it directly included gating and selectivity into the 
model. Lopreore et al. (2008) [13] developed a finite-volume-based technique to solve PNP 
in three dimensions, which decomposes the domain using a dual Delaunay-Voronoi mesh. 
Neuen (2010) [16] developed a semi-implicit finite element-based scheme to simulate three- 
dimensional, multi-scale extended PNP. Gardner and Jones (2011) [6] simulated a potassium 
channel modelled with PNP in two dimensions using a finite difference method with TR- 
BDF2 time integration. Much of the numerical schemes in [B] is based on the previous 
work [7], a one-dimensional model of the same channel. Hyon et al. (2011) [TT] presented 
another finite element method with back-Euler method in time to investigate the effects of 
finite size of the ions by modifying the PNP via introducing a repulsive potential energy into 
the total energy. Horng et al. (2012) [10] applied the multiblock Chebyshev pseudospectral 
method and the method of lines to solve a one-dimensional modified PNP modelling the 
finite-sizeness of the ions via a local model. 

One of the characteristics of the nonlinear PNP equations is that its overall behavior 
is very sensitive to the boundary conditions. [9] This presents a challenge for accurate and 
efficient numerical simulations, as generally the boundary conditions will have to be dis- 
cretized and approximated. In this paper, we shall investigate the effects of discretization 
error on the Poisson-Nernst-Planck equations, in particular discretization of the boundary 
conditions and the equations at the boundaries. We will demonstrate that the conservation 
properties of the numerical methods could be critical in obtaining the long-time behavior of 
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the solutions. 

The paper is organized as follows. We start by defining and simplifying the equations we 
are working with, in Sec. |2} including the introduction of the quantities that shall be preserved 
by our numerical schemes: the total concentration of each ion species in Sec. 12.11 and the 
energy dissipation law in Sec. 12.21 We then describe our numerical schemes in Sec. El which 
presents an approach to conserve the total ion concentrations exactly and approximate the 
energy dissipation law closely. Finally, we shall discuss the results of simulating the system 
using our numerical schemes in Sec. HJ 



2 Governing Equations 

Consider the PNP equations [3 [7] 
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,i = l,2,...,N, 
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(2) 



where q is the ion density for the i-th species, is the diffusion constant, Zi is the valence, 
e is the unit charge, ks is the Boltzmann constant, T is the absolute temperature, e is the 
permittivity, <j) is the electrostatic potential, po is the permanent (fixed) charge density of 
the system, and N is the number of ion species. [H] The equations are valid in a bounded 
domain Q with boundary dQ and for time t > 0. 

In this work, we shall use the no-flux boundary condition for Eq. (pQ). This may correspond 
to modelling the interior conditions of a channel that is in an occluded state, with closed 
gates at either end. Simulations of channels such as the KirBacl.l channel in such a state 
have been conducted in the past [3J. We shall use the Robin boundary condition for the 
Poisson equation, which models the effects of making the source of the potential across the 
channel partially removed from the ends of the channel. The formula for the boundary 
conditions are 

Vdd> -n = 0, i = 1,2,..., N, (3a) 
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for points on the boundary x G dQ. 

For some situations, such as a generic potassium channel separating potassium and 
chloride ion baths, the experimental data can be well-approximated by a one- dimensional 
model. [7j In one dimension, the equations (CEJ) and (J2J) are simplified as 
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for — L < x < L and t > 0, where L is the half of the length of the ion channel. The 
corresponding boundary conditions are 
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dx ksT 1 dx 



0. 



i ± )±r/-- = 0, forx = -L,L. (6) 

dx 



2.1 Total Concentration 

The total concentration per ion species is given by 

■■L 



i,tat(t) 



Ci(x, t) dx, z = l,2,. 
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(7) 



Due to the no-flux boundary conditions ()6]), the total concentration of each ion species 
is constant in time. This can be verified easily by differentiating (0) with respect to time, 
then applying the convection-diffusion equation (j3J) and no flux boundary condition OH]). 

One of the metrics we can use to evaluate different numerical schemes is therefore to 
measure how well the total concentration is conserved in numerical simulation. Ensuring 
that the total concentration for each species Cijot is constant will be the idea behind the 
schemes presented in this work. As will be seen in Sec. IU the preservation of the conservation 
property is crucial for producing correct numerical results over long time scales. 



2.2 Energy Dissipation 

The governing equations (jlj) and §5§ for the transport of ions can be derived from the energy 
of the system using variational principles. Similar to [TTJ, the total energy for our specific 
system is defined by 
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dx + -(0 + 0(L) + 0„0(-L)), (8) 
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where c^o are constants called "reference concentrations". Using the Poisson equation ([5]) 
the total energy can be written as 
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(9) 
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where the last term is the contribution of the electric energy from the boundaries. The total 
energy E satisfies the energy dissipation property 
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where /ij is the chemical potential of i'th ion species defined by the variational derivative of 
the energy with respect to the concentration Cj 
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The energy dissipation law f llOp can be derived by taking the time derivative of the total 
energy (|S]) and applying integration by parts, Eqs. (jl])-(JSD an d the boundary condition (J5]): 
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The rate of energy decay ( flOl) can be obtained by using the boundary condition (|6]) to show 
the last two terms on the RHS of ffT2l) cancel each other. 



2.3 Parameters and Nondimensionalization 

We specify the units and the parameters using the approximate values corresponding to the 
KcsA potassium channeljl]. In our ID model, the cylindrical channel takes a diameter of 10 A 
and a length of 120 A. We shall assume no permanent charges or selectivity for the purposes 
of this simulation. We consider the case of two ion species, i.e. N = 2, with the initial 
concentration for each ion being 2 molar, resulting in an initial number density (number of 
ions per unit volume) of 1.2044 x 10~ 3 ions/A 3 . The combination of the parameters ksT/e is 
approximately 0.025 V, assuming the temperature is T = 298 K. The permittivity e = e r eo 
is determined by the value of the vacuum e = 8.854187817 x 10~ 12 F/m and the relative 
permittivity e r (78.5 for water). 

The values of the diffusion coefficients Di depend on both the ion species and the channel. 
The only net effect of different diffusion constants is the rate of evolution of the system. 
Typical values for the diffusion coefficients for ion species in a channel are around 10 9 A 2 /s.[8] 
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We will select both diffusion coefficients to be equal to each other, causing them to take a 
value of one after nondimensionalization. 

The parameter 77, as a component of the Robin boundary condition (I3bl) . is an aggregate 
of multiple physical constants and is highly dependent on the properties of the surrounding 
membrane. Modelling the experimental setup as an electrical circuit shows that the quantity 
Aei/rj, where A is the area of the membrane and q is the permittivity of the membrane, has 
units of capacitance and is related to charge storage. The most significant charge storage 
contributing to Aei/rj is in fact the membrane capacitance, so we may surmise that the 
primary contributor to 77 is the membrane capacitance. If a very high capacitance to ground 
is present, 77 is approximated by the appealing formula 77 = Aei/C, where C is the capacitance 
of the membrane, however realistically 77 is much smaller than that. In this work, we shall 
take 77 = 2.78 x 10~ 3 A for our numerical simulations, but will also examine the effects of 77 
over a range from 10~ 5 A to 60 A. Changing the value of 77 might correspond to adding a 
parallel capacitance in experiment. 

Define the dimensionless variables and parameters c- = q/co, x' = x/L, t' = t/(L 2 /D ), 
D'i = Di/Do, <f)' = <ft/4>o, where Co is the average of the initial charge concentration, L is 
the half of the channel length or computational domain, Dq is a typical diffusion coefficient, 
0o is a characteristic value of the electrostatic potential such as the boundary value. Then, 
non-dimensionalizing the Nernst-Planck Eq. (j4j), we obtain 



df dx' ' 1 
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where Xi '■= ^o/^bT. 
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From the above, the dimensionless parameter xi 
ized Poisson Eq. (JSJ) is given by 
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Here, the dimensionless parameter e' is defined as e' := e/et where et is the characteristic 
permittivity chosen to be the value for water: e t = 6.950537436 x 10 -20 F/A. The non- 
dimensional parameter X2 is approximately 125.4 with these values. The corresponding 
dimensionless boundary conditions are 



D' 



dx' +Xl r 'dx' 
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+ 77' 



dn 



0, for x = —1, 1, 



(15) 



where 77' := fj/L. 

We drop the primes when we present our numerical methods for clarity. 

3 Numerical Methods 



We present a method for deriving numerical schemes that would conserve total concentra- 
tion of each ion species exactly if computations were performed without round-off errors. 
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We will illustrate the method by describing a mass-conservative scheme (i.e. preserving ion 
concentration exactly) for solving the nonlinear systems of PDEs (fT3l) and (I14p . The ex- 
tension of the method to the multi-dimensional case is straightforward. This scheme uses 
the trapezoidal rule and the second-order backward differentiation formula (TR-BDF2) in 
time and the second-order central differencing in space. The TR-BDF2 scheme is implicit 
in time, resulting in a system of nonlinear equations after discretization. Instead of using 
the Newton-Raphson method for solving the large nonlinear systems at each time step, we 
present a simple iterative scheme which is easy to implement and can solve the systems 
efficiently 



3.1 Discretization in Time 

For time-stepping, we shall use a slight modification of the scheme described in [TJ, which 
combines the trapezoidal rule with the second-order backward differentiation formula. 

(1) TR step: 

cr +7 ' fc+1 -7^/(cr 7 ' fc+1 ,0" + ^) =c? + 7 ^/(cr,0«), ^ = 1,2,^ = 0,1,2,..., 

d_ ( a<ft"+7.fc+l S _ _ / poL2 2 n+li k+l\ (16) 

dx dx J ~ \ <j> e t ~T A2 Z^i=l 

(2) BDF2 step: 

„w+l,i+l 1-7 A* f( r n + 1 > l + 1 A,n+l,l\ 1 n+7 (1~7) 2 r n ,* _ 1 9 7 — 1 9 

-i 2-7 "J i i / ~~ 7(2-7) i 7(2-7) « ' — > > — ' > > ' ' ' ' 

ax ^ fc ax y — ^ oet r a2 Z-(i=i J> 

(17) 



where /(q, 0) is defined as the right-hand side of (IT5|1 

= £ {a 

We take 7 = 2 — v2, which minimizes the local truncation error. j7] 

Removing the inner iterations, corresponding to the indices k in (fTB"|) and / in (JT7J), 
Eqs. ffT6"|) and ([17]) is the TR-BDF2 scheme requiring a nonlinear solver for the two systems 
of nonlinear equations: (1161) for (c n+7 , 0™ +7 ) at the grid points and (I17p for (c n+1 ,0 n+1 ). 
With the inner iterations, Eqs. (fTB"]) and (TT7j) provide a simple iterative scheme for solving 
the systems of nonlinear equations. For instance, at k-th iteration, we update the array 
c n+7,fc+i a ^ ^he g r j^ p i n ts by solving the first equation of (1161) which is a tri-diagonal system 
after the spatial discretization, since the values of <fi n+ ^> k are known at fc-th iteration; then, 
we update n+ T' fc+1 using the second equation of (fTBl) . We perform the inner iterations until 
convergence and, as shown later, choosing two inner iterations k = 2 and I = 2 would be 
sufficient. As for initial guesses at the n-th time step, we choose ^ n+7 '° = <fi n for fflBT) and 
0n+i,o _ 0n+7,/«+i £ Qr w i^h k corresponding to the last inner iteration at the previous 
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inner iteration. As shall be seen in Sec. HJ without any such inner iterations (k = I = 0), 
one could only attain first-order accuracy in time; on the other hand, with just one inner 
iteration (k = I = 1), one can attain second-order accuracy in time. In other words, the 
simple iterative scheme is very effective in solving the systems of nonlinear equations. 



3.2 Discretization in Space 

Next, we provide the discrete equations for the spatial differential operators in Eqs. ( fl6l) 
and ( I17jl . Let's divide the dimensionless interval [—1, 1] to J subintervals, Xj = —1 + jAx, 
where Ax = 2/ J and j — 0, 1, • • • , J. We denote the numerical values of g(x, t) at (xj, t n ) 
by and g(x) at Xj by gj. We present the standard second-order central differencing 
schemes for the spatial differential operators here to facilitate the description of the mass- 
conservative scheme which depends on the details of the discretization at the interior grid 
points (-J + 1 < j < J - 1). 

The ion diffusion term in Eq. ( [13]) is discretized as 

TA Di ^) [x,) m> • (19) 

The term driven by the electrostatic potential gradient in Eq. (fT5|) is given by 

d ( D fy\ t x \ w A,j+iQ, i+ iO i+2 - <pj) - Ih.j \<:.j \{Oj - (t>j,j- 2 ) ^ 
dx V 1 dx ) 4(Ax) 2 



The Laplacian in the Poisson Eq. (TH|) is approximated by 
d 



dx \ dx J ] (Ax) 



e i+ i0 i+ i - (e j+ i + e 3 _i)<Pj + e^i^-i • (21) 



3.3 Discretization of Boundary Condition 

We shall implement the boundary conditions using two different schemes. The first scheme is 
obtained by applying standard finite differencing to the boundary conditions, and the second 
is obtained by requiring the conservation of ions within the channel. As shown later, it is 
critical to preserve the ion concentrations for accurate numerical solutions. 



Standard Implementation 

Applying the forward differencing to the right-hand side of the Nernst-Planck equation ffl~3l) 
at the left boundary and using the no- flux boundary condition in ffl5l) . we obtain 
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Di > 1 %AxY (22) 
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It is similar at the right boundary. We implement the Robin boundary condition in (|T5|) 
with the second-order central differencing using ghost grid points as 

(0o -0-) - Z 1 ,/'' = 0, implying 0_i = 0! - ^^(0o - 0-), (23) 
2Ax 7] 

and similarly <p J+1 = $j_ x ((f) j - <f> + ). 

V 

Conservative Scheme: TR Step 

The no-ffux boundary condition in (fT5|) implies that the total concentration of each ion 
species is constant throughout time. Thus, we discretize the equations by requiring the 
numerical value of the total concentration be conserved exactly in time. 

First, we approximate the total concentration Cijot{t n ) defined in Eq. ([7]) using the trape- 
zoidal rule as follows 

c ttot = E c h Ax + ( c "o + ch) (24) 
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Let us examine the change of the total concentration in the TR step f)16p . 
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- t ,J i,J (25) 

7At 7At / v ; 

This summation has a telescoping effect where most of the interior terms cancel each other 
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and we are left with 



n+7 _ n 
u i,tot ^i,tot 
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We can achieve the conservation of the total concentration c™^J = c" tot , if we discretize the 



Nerst-Planck equation ([131) at the left boundary 
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and at the right boundary 
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It is important to point out that Eq. ( |2"7|) can be seen as discretizing Eq. ffl3l) using a 
first-order finite difference with grid size Ax/2 and using the no-flux boundary condition 
(TT5|) . Eq. fl2T|) can be rewritten as 
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Conservative Scheme: BDF2 step 

We can rewrite Eq. ffT7|) in such a way that the numerical value of the derivative of the total 
concentration becomes a linear combination of the result from the TR step and the right 
hand side of equation (fl~3j) evaluated at the n + 1th time step. 
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n+l _ n+7 1 _ n+7 _ n 

(1-7) At "2-7 7 At + 2-7 /lj j [ ' 

As with the TR step, almost all of the interior terms cancel in a telescoping sum, and we 
can require the exact conservation of the total concentration c"^,J = c™^,* in order to obtain 
the discretization of the Nernst-Planck equation (II 3D at the boundaries for the BDF2 step: 



n+ l n+7 1 /n+7 „n \ H x (V^ 1 _ r n + 1 ^ 

c i,o c i,o 1-7 / c i,o ~ c i,o \ 2 ^i±KH,i c i,o > 



;i - 7) At 2 - 7 V 7 At J 2-7 (Ax) 2 

Xi* A,ocy (g - gl) + A,!^ 1 (g - 0g) 
+ 2- 7 2(Ax) 2 ' Ki) 



n+l _ n+7 



7 )At 




7 2(Ax) 2 



n 7 1) 

(32) 



Equation f[3"Tj) can be seen as discretizing only the term /(c"^ 1 ) in Eq. (fT7|) using forward 
difference with grid size Ax/ 2 and using the no- flux boundary condition in 015p . Eq. 
can be viewed similarly at the right boundary. 



4 Numerical Results 

4.1 Validation and Convergence Results 

To validate the accuracy our numerical method, we compare the steady-state solution from 
our dynamic simulations of PNP with that of the Poisson-Boltzmann solution taken from the 
work [12] . Figure [1] shows that our steady-state solutions match perfectly with those in [12] 
for two sets of parameters: one with rj = e = 2~ 2 and the other 77 = e = 2~ 6 while keeping 
the other parameters constant: </>_ = — 1, (f) + — 1, D\ — D2 — 1, Xi — 1> X2 = 5") and po = 0. 
The maximum difference in between the two solutions is less than 5.6 x 10~ 5 . To get the 
steady-state solution, we have used the mass- conservative TR-BDF2 method described in 
previous sections with 2048 grid points in the interval [—1,1] as in [12] and the time-step 
size 10 -4 . At time t — 0, the initial profiles for the ion concentrations are uniform in space. 
In this case, our time- dependent solution is close to the steady-state solution for the time 
t > 2. We have also verified that our solutions agree with those in [12] for other sets of 
parameters as well, although they are not shown here. 

We have also checked the orders of convergence of our methods. The discretization 
method described in the previous section always has 0(Ax 2 ) convergence in space, regardless 
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Steady State Comparison 




Figure 1: Comparing our steady-state solution (the dashed lines) using TR-BDF2 method 
with that of the Poisson-Boltzmann equation (the solid lines) obtained in [12]. The param- 
eters are e = 2~ 2 , 2~ 6 , 77 = e, 0_ = —1, 0+ = 1. 



whether we have implemented the mass-conservative difference scheme or not. The order of 

, . , r , , I $(2 As) - $(4Ax)| , , . . 

convergence m space is computed using the formula log 2 - — — — — - — — , where <P( Ax) 

" |<P(Ax) — <P(2Ax)| 

denotes the numerical solution of the potential at the point (x, t) = (0.904,0.02) obtained 
with the spatial resolution Ax. In this case, the time step size is chosen to be very small 
At = 10~ 6 so that the discretization error is dominated by that in space. 

To obtain the numerical orders of convergence in time, we compute the numerical solu- 
tions with three different time-step sizes At, 2 At and 4At and then calculate the numerical 
order of convergence p by computing the ratio $(2At) — $(4At))/($(At) — $(2At)) at the 
fixed position and time (x,t) = (0.904,0.02). Here, the spatial resolutions in these simula- 
tions are kept the same, Ax = 0.002. The numerical convergence results in time are given 
in Table [TJ We find that, if one did not perform inner iterations (k = in ( fl6|) and / = in 
( TT7|) ). the convergence of TR-BDF2 would be first-order in time. If we include at least one 
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inner iteration (fc > 1 and I > 1), then the convergence becomes second-order as expected. 



At 


5 x kt 5 


2.5 x 10~ 5 


1.25 x 10" 5 


order of convergence for TR-BDF2, no inner loops 
order of convergence for TR-BDF2, two inner loops 


1.0016 
2.2197 


1.0008 
2.1779 


1.0028 
2.2143 



Table 1: The numerical order of convergence in time for the mass-conservative TR-BDF2 
method solving the PNP equations in one dimension for two ion species. The non- 
dimensionalized physical parameters are e = 1, rj = 4.63 x 10~ 5 , 0_ = 1, <p + = —1. The 
calculations are performed with Ax = 0.002 and the numerical solution of 4> is evaluated at 
the point (x,t) = (0.904,0.02). 



4.2 Evolution of the Distributions of the Ions 

First, we examine the evolution of the ion concentrations and the electrostatic potential 
starting from a uniform ion distribution of two ion species of opposite valence Z\ — 1 and 
z 2 = —1: Ci(x, 0) = 1, i = 1,2, for — 1 < x < 1. The prescribed electrostatic potentials 
on the left and the right at far-field are 0_ = 1 and <p + = — 1 respectively. The physical 
parameters are specified as in Sec. 12.31 In the rest of this work, unless we specify otherwise, 
the non-dimensionalized parameters are chosen as D\ = D 2 = l,Xi = 3.1, x 2 = 125.4 and 
t] = 4.63 x 10~ 5 , as they were defined in Sec 12.31 Due to the symmetries of the initial and 
boundary conditions, the parameters and the domain, the profiles for the concentrations of 
the two ion species at any time are symmetric with respect to the center of the channel, 
x = 0. 

Figure |2] shows the profiles of the ion concentration with the valence z 2 = — 1 and the 
electrostatic potential at the times t = 0, 0.01, 0.05, and 1. The Robin boundary condition 
( TL5"]) for the electrostatic potential drives the ions with negative charges toward the left 
boundary and the no-flux boundary condition (j!5p for the ions causes those charges to 
accumulate at the boundary. In this case, the ion concentrations keep their uniform profile 
in the bulk of the domain away from the two ends, while the electrostatic potential changes 
from an initially linear profile to one that is essentially constant (zero) except for the sharp 
gradient at each end. We find that the existence of the thin boundary layers requires high 
spatial resolution or small Ax in the simulation. The numerical results would be far away 
from the correct solution if we chose Ax > 0.05. These results show the overall behavior of 
the system as time elapses. 

4.3 Comparison between Mass-conservative and Standard Schemes 

Next, let us compare the numerical results from a standard discretization (called as the non- 
conservative schemes) of the boundary conditions, ( I22p . with those obtained from the mass- 
conservative schemes ( 12T|) and ( 131]) . Figure [3] shows the ion concentration profiles and the 
electrostatic potential at time t = 1 obtained from both the mass-conservative schemes(the 
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solid lines) and the non- conservative schemes (the dashed lines). The parameters in the 
computations are the same as described in the previous Sec. 14.21 To make fair comparison, all 
other aspects are kept same, including the time-step scheme (TR-BDF2), the discretization 
scheme for interior points of the domain, the initial condition, the physical parameters, the 
time-step size At and the space resolution Ax. As shown in Fig. [3(a), the ion concentration 
from the non-conservative scheme is substantially lower than that from the mass-conservative 
scheme and the variations near the boundaries are much smaller in the result from the 
non-conservative scheme. Furthermore, the electrostatic potential obtained from the non- 
conservative scheme, shown in Fig. [3(b), has a linear profile with non-zero slope in the middle 
of the domain and much milder slopes at the boundaries, when compared with that from 
the mass-conservative schemes. 

Because of the no- flux boundary conditions (I3al) . the total concentration of each ion 
species should be invariant in time. Figure H] shows that the mass-conservative scheme 
preserves the conservation of the ions perfectly (up to the level of roundoff error) over a 
long period of time, while the total number of ions at the time t = 1 obtained from the 
non-conservative scheme is reduced to less than half of the original amount. 

Figure [5(a) shows that the total energy E as a function of time t for both the conservative 

and non-conservative schemes. The total energy obtained from the mass-conservative scheme 

approaches the minimum energy state much faster than that from the non-conservative 

scheme. More importantly, in Sec. 12.21 it is shown that the total energy of the system E 

defined as (IE]) satisfies the energy dissipation law ( 1101) . In Fig. |5jb) , we plot the rate of change 

in energy, for the mass-conservative (the solid line) and the non- conservative schemes (the 

dotted line) obtained by using a second-order finite difference based on the numerical result 

E(t) shown in Figure [5(a). In the same graph, we also plot the expected dissipation rate 

given by the right-hand side of (11 01) . computed using the second-order central differencing 

and trapezoidal rule and shown by the dashed line for the conservative scheme and the 

dash-dotted line for the non- conservative scheme in Fig. [5(b). It shows that the numerical 

result from the conservative scheme (the solid line) agrees with the energy dissipation law (the 

dashed line) very well. In contrast, the corresponding results for the non- conservative scheme 

show that the energy dissipation law is not satisfied after a short period of time. This is due 

to the fact that the total concentration from the non-conservative scheme displays very poor 

performance in conserving the total concentrations. The results show that the discretization 

of the boundary conditions have profound impact on satisfying the physical properties: the 

energy dissipation law and the conservation of the total number of ions. 

In addition to energy decay, we compute the maximum rate of change in the concen- 
tre ■ 

trations of the species over the domain, i.e. max I — — — i . It is notable from the time 

i,-i<x<i dt 

derivative of concentration shown in Fig. [6] that the numerical results from the conserva- 
tive numerical scheme steadily approach the equilibrium in time. On the other hand, the 
non-conservative scheme is approaching a steady state much faster initially, but, later in 
time, the non-conservative scheme's behavior changes and it does not appear to reach a 
steady state. This result emphasizes the necessity of the conservative numerical scheme for 
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long-time simulation. 

4.4 Effect of Parameters 

The size of the difference in the results from conservative and non-conservative schemes 

ec L 2 

depends on the non-dimensional parameter \2 — —, • For the physical model of the 

Wi 

ion transportations, the value of X2 can be arbitrarily large, depending on the values of 
average ion concentration Co and the applied electrostatic potential 0o at the boundaries. 
Consequently, it is important to pay attention to the size of the dimensionless parameter \ 2 . 
In Fig. [31 we have shown that, for \ 2 = 125.4, the results of non-conservative schemes are 
far away from the correct results. Figure [7(a) and (c) show the profiles of the electrostatic 
potential at a fixed time t — 1 from both the conservative and the non-conservative schemes 
with two more different values of \2 = 31.35 and 501.6, while keeping all other parameters 
the same as those for Fig. [3] At t — 1, the system has reached the steady state, shown 
by the constant values for the conservative scheme in the energy plots of Fig. [7(b) and (d). 
Comparing the graphs of potential in Fig. [3(a), (c) and Fig. [7J we find that the value of X2 
primarily affects the width of the boundary layer, with larger \2 resulting in thinner boundary 
layers. A thinner boundary layer transitions much more sharply near the boundaries, and 
thus requires more computational grid points in the region and more truthful discretization 
of the boundary conditions. This causes the differences in electrostatic potential profiles and 
the energy dissipation in time (shown by Figs. [7(b) and (d)) between the conservative and 
non-conservative schemes to be greater as one increases \2- A thinner boundary layer also 
affects performance with regard to the energy dissipation law, which is not shown here in 
plots. Larger \2 leads to a larger discrepancy between the decay rate of the total energy 
(the left-hand side of Eq. [TUj) and the energy dissipation rate (the right-hand side of the law 
Eq. [TO]) . and this discrepancy gets worse faster for the non-conservative scheme than for the 
conservative scheme. 

Finally, we examine the effect of the parameter 77 in the Robin boundary condition (13bj) 
on the numerical results. As noted in Sec. 14.11 the steady state changes dramatically if the 
relative values of the physical parameters rj and e are changed. In order to determine the 
effect of rj itself on the results, we have tested a range of non-dimensionalized values for 77 
ranging from 10~ 6 to 1, while holding e at its constant non-dimensionalized value of 1. We 
find that, when 77 increases from 10 -6 to 0.001, the concentration profiles at the steady state 
do not change much, having a maximum relative difference of only 10~ 4 , but this property 
does not generalized to larger 77. We also find that the discretization error, especially for 
the non-conservative scheme, is significantly affected by the value of 77. For large values 
of 77, say 77 > 0.1, the growth of the discretization error of the non-conservative scheme is 
rather slow, and consequently the concentration and electric potential profiles obtained from 
the non- conservative scheme are close to those obtained by the mass- conservative schemes. 
An example of this property is shown in Fig. [HI It appears that, for 77 = 0.5, the total 
energy from the non-conservative scheme decreases linearly in time after an initial sharp 
drop, becoming negative at later time. On the other hand, the conservative scheme reaches 
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a steady state very quickly and does not deviate from it. For small values of r\ such as those 
shown in Fig. [3j both the conservation property of the total concentrations and the energy 
dissipation law deteriorate at a fast pace for the non-conservative scheme, and the difference 
between the results from the conservative and the non- conservative schemes grows bigger as 
77 gets smaller. 

5 Conclusion 

The primary objective of this work is to investigate the effects of conservation property of 
discretization schemes on the numerical results. We have shown that, with regard to the 
PNP equations, whether a numerical method preserves the mass conservation could have a 
critical impact on the behavior of the system, especially the steady state results. We have 
provided a discretization scheme that preserves the mass conservation exactly (excluding the 
round-off errors) and the energy dissipation law well for long-time simulation. 

Our method is implicit in time and second-order accurate in both space and time. We 
have verified that approximating the fully implicit solution is necessary for second-order 
convergence in time. Further, we find that one can avoid using Newton-type nonlinear 
solvers by performing a simple iterative scheme. 

In this work, we have simulated the equations with realistic physical parameters, partic- 
ularly investigating the effect of the non-dimensional parameters \2 in the Poisson equation 
and 77 in the Robin boundary condition for the electrostatic potential. We find that the 
mass- conserving scheme is more robust to changes in parameters, especially changes to the 
value of 77. 

Although this work makes good progress in constructing an accurate method for solv- 
ing the Poisson-Nernst-Planck equations numerically, there are many challenges remaining. 
First, one of them is to account for the finite size of the ions as its effect is enormous 
considering the narrow width of the ion channels. |1 11 [TO] Second, for most ion channels, 
the appropriate boundary conditions are Dirichlet-type. We will investigate the possibility 
to preserve the energy dissipation law exactly instead of the mass and study the effect of 
the conservation on long-term behavior of the simulation. Third, we would like to include 
distributions of permanent charges for studying selectivity of ion channels. 
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(b) Electric Potential 
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Figure 2: Simulation results using the mass- conservative TR-BDF2 method for e = 1, rj = 
4.63 x 10 -5 , 0_ = 1, cf) + = — 1. The calculations were performed with At = 10~ 4 and 
Ax = 0.002. (a) The concentration profiles for the ion species with the valence Z2 = — 1, 
C2(x,t), are plotted at the times t — (the solid line), 0.01 (dashed), 0.05 (dotted) and 
1 (dash-dotted), (b) The corresponding time sequence of the electrostatic potential <fi is 
plotted. 
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(b) Electric Potential 
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Figure 3: Comparison between the simulation results from the mass-conservative and the 
non-conservative schemes for e — 1, rj — 4.63 x 10~ 5 , 0_ = 1, <ft + = —1, T — 1. The 
calculations were performed with At = 10~ 4 and Ax = 0.002. (a) The ion concentration 
profiles of c 2 from the mass-conservative method (the solid line) and the non-conservative 
method (the dashed line), (b) The corresponding electrostatic potentials. 
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(a) Total Ions, Species 2 
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Figure 4: (a) The total ion concentration for species 2 as a function of time from the 
simulations using the mass-conservative (solid) and non-conservative (dashed) schemes, (b) 
The relative error in total concentration for both species. The parameters are identical to 
those in Fig. |3j 
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(b) Energy Law 
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Figure 5: (a) The total energy as a function of time from the simulations using the mass- 
conservative (solid) and non-conservative (dashed) schemes, (b) The rate of change in energy, 
obtained from the graph (a) and the right-hand side of Eq. ( jTUl) . The solid and the 
dotted lines correspond to the left-hand side of Eq. (jTUj) for the mass-conservative and the 
non-conservative schemes respectively. The dashed and the dash-dotted lines correspond to 
the right-hand side of Eq. (fTOj) for the mass- conservative and the non- conservative schemes 
respectively. The parameters are identical to those in Fig. [3J 
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Figure 6: The maximum rate of change in ion concentrations as a function of time for the 
non-conservative(the dashed line) and conservative (the solid line) schemes. The parameters 
are identical to those in Fig. [3j 
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(a) Electric Potential, % =31.3499 
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(c) Electric Potential, % =501.5992 
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(d) Energy Total, % =501 .5992 
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Figure 7: Comparison between the simulation results from the mass-conservative and the 
non-conservative schemes for different values of the non-dimensional parameter xi- The 
calculations were performed with At = 10 -4 and Ax = 0.001. The other parameters are 
identical to those in Fig. [3] 
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(b) Electric Potential 
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Figure 8: Comparison between the simulation results from the mass-conservative and the 

non-conservative schemes for 77 = 0.5. The other parameters are identical to those in Fig. [3j 

(a) The ion concentration at the non-dimensionalized time T = 1. (b) The electric potential 

at the non-dimensionalized time T — 1. (c) The change of the total energy in time. 
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